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APPLICATIONS OF THE LATTICE BOLTZMANN METHOD TO COMPLEX AND 

TURBULENT FLOWS 

LI-SHI LUO*, DEWEI Qlt, AND LIAN-PING WANG* 

Abstract. We briefly review the method of the lattice Boltzmann equation (LBE). We show the three- 
dimensional LBE simulation results for a non-spherical particle in Couette flow and 16 particles in sedimen- 
tation in fluid. We compare the LBE simulation of the three-dimensional homogeneous isotropic turbulence 
flow in a periodic cubic box of the size 128 3 with the pseudo-spectral simulation, and find that the two 
results agree well with each other but the LBE method is more dissipative than the pseudo-spectral method 
in small scales, as expected. 

Key words, lattice Boltzmann method, turbulent flow, 3D homogeneous isotropic turbulence, spectral 
method, non-spherical particulate suspensions 

Subject classification. Fluid Mechanics 

1. Introduction. More than a decade ago, the lattice-gas automata (LG A) [5, 24, 6] and the lattice 
Boltzmann equation (LBE) [17, 12, 2, 22] were proposed as alternatives for computational fluid dynamics 
(CFD). Since their inception, the lattice-gas and lattice Boltzmann methods have attracted much interest 
in the physics community. However, it was only very recently that the LGA and LBE methods started 
to gain the attention from CFD community. The lattice-gas and lattice Boltzmann methods have been 
particularly successful in simulations of fluid flow applications involving complicated boundaries or/and 
complex fluids, such as turbulent external flow over complicated structures, the Rayleigh-Taylor instability 
between two fluids, multi-component fluids through porous media, viscoelastic fluids, free boundaries in flow 
systems, particulate suspensions in fluid, chemical reactive flows and combustions, magnetohydrodynamics, 
crystallization, and other complex systems (see recent reviews [3, 16] and references therein). 

Historically, models of the lattice Boltzmann equation evolved from the lattice-gas automata [5, 24, 6]. 
Recently, it has been shown that the LBE is a special discretized form of the continuous Boltzmann equation 
[8, 9]. For the sake of simplicity without loss of generality, we shall demonstrate an a priori derivation of 
the lattice Boltzmann equation from the continuous Boltzmann equation with the single relaxation time 
(Bhatnagar-Gross-Krook) approximation [1], The Boltzmann BGK equation can be written in the form of 
an ordinary differential equation: 

where D t = dt + , f = f(x, t) is the single particle distribution function, A is the relaxation time, 

and is the Boltzmann distribution function in ^-dimensions, in which p, u and 9 = k/jT /rn are the 
macroscopic density of mass, the velocity, and the normalized temperature, respectively, T, ks and m are 
temperature, the Boltzmann constant, and particle mass. The macroscopic variables are the moments of the 

*ICASE, Mail Stop 132C, NASA Langley Research Center, 3 West Reid Street, Building 1152, Hampton, Virginia 23681- 
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Virginia 23681-2199. 
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distribution function / with respect to the molecular velocity 


P= J f<% = f f m d£, (1.2a) 

P* = jtfdt = jtfWdi, (1.2b) 

pO=\j{£-u) 2 fdZ=\j{ti- uf /<») da . (1.2c) 

Equation (1.1) can be formally integrated over a time interval Sp. 

f(x + as t , a, t+s t ) = e~ St,x fix, a, t) + \ e ~ St,x J e</A / (0) (* + 4*\ 4, * + *') dt 1 . (1.3) 

Assuming that is small enough and /(°) is smooth enough locally, and neglecting the terms of the order 
0(( 5|) or smaller in the Taylor expansion of the right hand side of (1.3), we obtain 

fix + aSt, a, t + S t )~ f(x , 4, f) = --[/(*, £, £) - /^ 0) (a;, £, £)] , (1.4) 

T 

where t = A/<S t is the dimensionless relaxation time. The equilibrium / ,0) can be expanded as a Taylor 
series in u up to u 2 


/ (eq) = 


(2tt 0)W2 exp l 28 


l + ^ + (£ • w ) 2 


21 


8 


28 2 


u 

28 


(1-5) 


To obtain the Navier-Stokes equations, the hydrodynamic moments (/>, pu , and p8) and their fluxes 
must be preserved in finite discretized momentum space {£ a |a = 1 , 2 ,..., 6}, i.e., 


C¥ cx 

(1.6a) 

p U = J2Saf a = Y,^fa q) ’ 

a a 

(1.6b) 

PO = \ - «) 2 fa = l'£tt°~ U ^ /a e<5) > 

(1.6c) 


a a 


where f a = f a (x, t ) = W a f(x, £ a , t ) [8, 9]. It turns out that these moments can be evaluated exactly in 
discretized momentum space by using Gaussian-type quadrature [8, 9, 23]. 

We can derive the nine-velocity athermal LBE model on a square lattice in two-dimensions 

f a {Xi +e a 8 t , t + S t ) - fa(Xi, t) = —[fa(Xi, t ) - f^ q> (Xi, t)] , (1.7) 

T 

where the equilibrium fa q \ the discrete velocity set {e a }, and the weight coefficients {w a } are given by 



+ 
i— i 

CE 

Q 

B 

3(e a • u) 
c 2 

9(e a • it) 2 3u 2 ) 

2c 4 2c 2 J 

* 1 



(1.8a) 


f (0,0), 


a = 0 , 


f 4/9, 

a = 0, 


€-a — ^ 

(±1, 0)c, 

(0, ±l)c, 

a = 1, 2, 3, 4, 

w a = < 

1/9, 

a = 1, 2, 3, 4, 

(1.8b) 


{ (±T ±1)a/2 c, 

a = 5, 6, 7, 8, 


l 1/36, 

a = 5, 6, 7, 8, 



and c = S x /St- Equation (1.7) involves only local calculations and uniform communications to the nearest 
neighbors. Therefore it is easy to implement and natural to massively parallel computers. 
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The (incompressible) Navier-Stokes equation derived from the above LBE model is: 

pdtu + pu\7u = —VP + pvV 2 u , (1.9) 

with the isothermal ideal gas equation of state, the viscosity, and the sound speed given by 



It should be noted that the factor —1/2 in the above formula for v accounts for the numerical viscosity due 
to the second order derivatives of f a . This correction in v formally makes the LG A and LBE methods second 
order accurate. Similarly, we can derive the six-velocity and seven-velocity models on a triangular lattice in 
two-dimensions, and the twenty-seven- velocity models on a cubic lattice in three-dimensions [9] . 

There have been some significant progress made recently to improve the lattice Boltzmann method: (i) 
the generalized lattice Boltzmann equation with multiple relaxation times which overcomes some shortcom- 
ings of the lattice BGK equation [14]; (ii) use of grid refinement [4] and body-fitted mesh [10, 7] with inter- 
polation/extrapolation techniques; (iii) adaptation of unstructured grid by using the finite element method 
or the characteristic Galerkin method; (iv) application of implicit scheme for steady state calculation and 
multi-grid technique to accelerate convergence (see a recent review [16] for further references). 

In what follows we shall demonstrate the applications of the LBE method to simulate the flow of non- 
spherical particulate suspensions in fluid and homogeneous isotropic turbulence in a periodic box. 

2. LBE Simulation of Flows of Non-Spherical Particulate Suspensions. The flow of particulate 
suspensions in fluid is difficult to quantify experimentally and to simulate numerically in some cases. Yet 
the flow of particulate suspensions is important to industrial applications such as fluidized beds. There have 
been some successful simulations of the flow of spherical suspensions by using conventional CFD methods, 
such as the finite element method. However, the simulation of the flow of non-spherical suspensions still 
remains as a challenge to the conventional CFD methods. Recently the LBE method has been successfully 
applied to simulate the flow of non-spherical suspensions in three-dimensions [19, 20]. The success of the 
LBE method to this problem relies on the fact that the LBE method can easily handle the particle-fluid 
interfaces [15], and accurately evaluate the force on the particle due to the fluid flow [18]. 

We first simulate a single non-spherical particle in the Couette flow. The equilibrium states in a non- 
spherical particulate suspension in a 3D Couette flow are simulated for a particle Reynolds number up to 
320. Particle geometries include prolate and oblate spheroids, cylinders and discs. We show that the inertial 
effect at any finite Reynolds number qualitatively changes the rotational motion of the suspension, contrary 
to Jeffery’s theory at zero Reynolds number [13]. At a. non-zero Reynolds number, a non-spherical particle 
reaches an equilibrium state in which its longest and shortest axes are aligned perpendicular and parallel 
to the vorticity vector of the flow, respectively. This equilibrium state is unique, dynamically stable, fully 
determined by the inertial effect, the maximum energy dissipation state. Systems of either fifty cylinders or 
fifty discs in Couette flow are also simulated. Multi-particle interactions significantly change the equilibrium 
orientation of solid particle. The effect is stronger for cylinders than for discs. The details of this work will 
be reported elsewhere [21]. 

Figure 1 shows a 3D LBE simulation of sixteen cylindrical particles falling under the influence of gravity. 
The left figure illustrates the time evolution of the entire system of sixteen particles, while the right figure 
demonstrates the formation of inverted T configurations in the sedimentation, which has been observed 
experimentally. To the best of our knowledge, this phenomenon was first reproduced numerically by the 
LBE direct numerical simulation [20]. 
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Fig. 1. 3D LBE simulation of particles sedimentation in fluid. Particle size is D = 12 and L = 24. System size is 
N x x N y x N z = 140 x 150 x 35. The averaged single-particle Re « 16.9. (left) Evolution of 16 particles (from left to right and 
top to bottom), (right) Formation of inverted T configurations which are also observed in experiment. 

Table 3.1 

Parameters in lattice Boltzmann and pseudo-spectral simulations: L is the length of box side; TV 3 is the system size; v is 
the viscosity; u 1 is the RMS fluctuation of the initial velocity field; dt is the time step size; T is total integration time, Re A is 
the Taylor microscale Reynolds number; and M is the Mach number. 


Method 

L 

N 3 

V 

u' 

dt 

T 

Re A 

M 

Spectral 

2tt 

128 3 

0.01189 

0.993311 

0.002 

2 

35.0 

0 

LBE 

128 

128 3 

0.009869 

0.040471 

1 

1000 

35.0 

0.0687 


3. LBE Simulation of 3D Homogeneous Isotropic Turbulence. Homogeneous isotropic turbu- 
lence in a three-dimensional periodic cubic box remains as a stand problem in the field of direct numerical 
simulation of turbulence. Due to the simplicity of the boundary conditions, the pseudo-spectral method can 
be easily used to simulate the flow. Because of its accuracy, the pseudo-spectral result is often used as a 
benchmark standard. Here the LBE simulation of the flow is compared with the pseudo-spectral simulation. 

The parameters of the simulation are given in Table 3.1. The initial condition is a random velocity field 
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Fig. 2. LBE vs. Pseudo-spectral DNS of 3D homogeneous isotropic turbulence. System size is 128 3 . Re* = 35. (a) The 
energy spectrum E(k ) as a function of time, (b) The decay of the mean kinetic energy K and dissipation rate e. The results 
from the LBE simulation are scaled according to the dimensions used in the spectral simulation. 


with a Gaussian distribution and a compact energy spectrum 


k 2 

E(k) oc — exp 

fco 


A 21 

ko 


The boundary conditions are periodic in three dimensions. The Taylor microscale Reynolds number is defined 


as 


Re* = 


2 K(t — 0) A u R msA 


where K(t = 0) = {«o/2)y = (3u 2 MS /2)y is the volume averaged kinetic energy (of the initial zero-mean 
Gaussian velocity field u 0 with RMS component u HMS ) , and A is the transverse Taylor microscopic scale: 


A = ^15vul MS /e, 


where e is the dissipation rate. 

Figure 2 shows the energy spectrum E(k) as function of time, and the time evolution of the mean kinetic 
energy K and dissipation rate e. The lattice Boltzmann results (symbols) are compared with the pseudo- 
spectral results (lines). The LBE results agree well with the pseudo-spectral results. Obviously the LBE 
method is more dissipative, especially at high wave numbers k > |fe m ax, where fc max = |JV, and N is the 
number of mesh nodes in each direction. This is because the LBE method is only second order accurate in 
space and time and thus more dissipative than the pseudo-spectral method. 

4. Conclusions and Discussion. The above simulations were performed on a Beowulf cluster of 
Pentium CPUs. For the simulation of the particulate suspension, the code consists two part: the lattice 
Boltzmann method for the fluid and molecular dynamics (MD) for the solid particles [19]. Even though 
the MD part of the code is not yet parallelized, the speed of the code still scales well with the number of 
CPUs up to 32 CPUs when the system size is 64 3 and with fifty particles. Presently we can easily simulate 
a system of a few hundred particles on our Beowulf system. 
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As for the simulation of the 3D homogeneous isotropic turbulence, the LBE code without optimization 
has the same speed as the spectral code with a Beowulf cluster of eight CPUs (about Is per time step). 
However, we do expect the LBE code will scale linearly with the number of CPUs, but not the spectral code. 

Our current research includes particulate suspension in fluid with high volume fraction of particles, vis- 
coelastic and non-Newtonian fluids, and forced or free-decay homogeneous isotropic turbulence in a periodic 
cube by using the lattice Boltzmann method on massively parallel computers. 
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